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MODELING  AND  SIMULATION  OF  AN  IMPLODING 
PLASMA  RADIATION  SOURCE 


I  INTRODUCTION 

The  use  of  an  imploding  gas-puff  plasma  as  a  plasma  radiation  source  (PRS)  is  of 
current  experimental  interest.  In  particular,  a  neon  PRS  has  been  used  as  a  load  for  the 
GAMBLE  II  puised-power  device  at  the  Naval  Research  Laboratory,  and  has  recently  been 
measured1  to  produce  approximately  2kJ  of  kilovolt  if -shell  X-rays  with  a  peak  load  cur¬ 
rent  of  1.2  MA.  In  addition  to  the  experimental  applications,  theoretical  interest  has  grown 
concerning  the  plasma  and  atomic  physics  involved  in  the  high-density  (~  1019cm~3)  pinch 
commonly  achieved  in  these  experiments:  the  plasma  instabilities  and  radiation  mechanisms 
are  among  the  complicated  processes  which  are  still  not  completely  understood.  Between 
the  experimental  and  theoretical  research,  there  is  a  need  for  simple  models  of  the  physics  of 
an  imploding  plasma  with  accompanying  radiation  that  can  be  implemented  in  a  tractable 
computer  simulation  of  not  only  just  the  PRS  itself,  but  of  the  entire  puised-power  device 
(generator,  transmission  line,  switches,  etc.)  as  well.  As  a  primarily  inductive  load  on 
the  remainder  of  the  circuit,  the  evolution  of  the  PRS  affects  to  some  extent  the  behavior 
and  efficiency  of  the  entire  system.  Such  a  system  simulation  is  important  for  investigating 
parameter  regimes  not  currently  available  to  the  experiment  in  order  to  determine  opti¬ 
mal  design  and  operating  parameters.  As  this  type  of  simulation  typically  requires  a  large 
number  of  computer  runs,  in  terms  of  speed  and  storage  requirements  the  use  of  sophisti¬ 
cated  one-  or  two-dimensional  hydrodynamic  models  of  the  plasma  coupled  with  complicated 

atomic  physics  algorithms  is  inappropriate.  Thus,  the  intent  of  the  model2,3  contained  in 
Manuscript  approved  March  27.  19S6. 


the  present  report  is  to  include  just  enough  physics  to  adequately  describe  the  evolution  of 
the  PRS  temperature-dependent  effects  such  as  radiation,  ionization  and  ‘‘bouncing”. 
etc.),  yet  computationally  efficient  in  order  to  simulate  the  operation  of  the  entire  device. 


In  Section  II,  we  will  state  our  assumptions  and  derive  the  equations  of  motion  which 
govern  the  imploding  plasma.  Section  III  discusses  the  details  of  the  procedure  used  in  the 
numerical  implementation  to  convert  part  of  the  kinetic  energy  gained  by  the  plasma  during 
the  run-in  phase  into  thermal  energy,  ionization  and  radiation  w’hen  the  plasma  collides  on 
axis.  The  atomic  physics  model  we  use  is  described  in  Section  IV.  The  implementation  and 
use  of  the  code  IMPLODE  is  detailed  in  Section  V,  along  with  numerical  results  compared 
with  experimental  observations  (from  GAMBLE  IT)  and  parameter  optimization  results. 
It  will  be  shown  in  Appendix  A  that  the  ordinary  differential  equations  evolved  by  the 
code  conserve  energy  in  that  the  power  delivered  to  the  PRS  load  by  the  external  circuit 
is  accounted  for  in  various  forms  of  stored  energy  and  radiation:  this  is  important  from  a 
design  energy  inventory  point  of  view. 


II  DERIVATION  OF  EQUATIONS 

The  geometry  of  our  m<  del  gas-puff  plasma  is  that  of  a  cylindrical  annulus  (Fig.  I)  of 
length  i  with  inner  radius  at  r  =  6  and  outer  radius  at  r  =  a  inside  a  cylindrical  cavity  of  wall 
radius  rw.  The  gas-puff  is  composed  of  a  plasma  with  ion  mass  m,  =  AmP  and  mass  /unit 
length  ft.  The  annulus  itseif  is  further  divided  into  two  concentric  annular  regions:  an  outer 
skin  of  skin  depth  6  and  a  larger  inner  core.  The  boundary  between  the  two  regions  is  at 
a  radius  r  =  s  (t.e„  x  =  a  -  6).  The  current  I  in  the  plasma  is  assumed  to  flow  only 
in  the  skin  of  conductivity  c.  and  the  current  density  J  =  Jz  is  assumed  to  be  uniform. 
.7  =  7/tr(o2  -i2).  The  two  regions  (core  and  skin',  are  allowed  to  have  separate  temperatures 
and  are  not  in  thermal  contact,  that  is.  there  is  no  thermal  diffusion  between  the  regions. 
Quantities  referring  to  the  skin  will  be  labelled  with  a  subscript  s  and  those  of  the  core  with 


J&.  '■  JL-<* 


a  subscript  a  thus,  the  temperature  and  volume  of  the  core  and  skin  are  respectively  Te,  Ve 
and  T..V.. 

We  begin  with  the  continuin’  equation 


-37  +v  •  Vp  =  -pV  •  v, 
ot 


where  p  is  the  mass  density, 


p  =  n,m,  +  n«me  =  n,mt  -f  Zn,me  =  nm,(l  +  Zme/mt)  «  nm,  =  — r— = — jtt,  (2) 

—  r ) 

and  the  number  density  n  is  defined  as  the  ion  density  n,.  The  first  assumption  will  be 
Al.  The  mass  density  (or  ion  density)  is  uniform  over  b  <  r  <  a  (i.e.,  the  same  in  both  skin 


and  core),  but  is  allowed  to  vary  in  time. 


Thus, 


n(t)  - 


m,-ir(a2(t)  —  b?(t)) 

Therefore,  with  Vn  =  0,  (1)  becomes  an  equation  for  the  velocity  profile  t»(r) 


n  ad  —  66  Id 

-  =  -2— — —  =  — 
n  a?  -  hr  .  r  dr 


Now,  the  second  assumption 

A2.  The  inner  radius  moves  with  the  same  speed  as  the  outer  radius. 


a  s  b  =  u. 


With  this,  the  equation  for  v(r)  is 


1  d  2u 

77b' 


which  has  the  solution  matching  the  boundary  condition  (5) 

abu  1  u  Ai 

v(r)  =  - - - - rT= - -  Ajf.  (') 

a  +  b  r  a  —  6  r 

Note  that  this  form  for  v (r)  also  holds  as  6  —  0,  in  which  case  A\  — *  0  so  that  »(r)  =  u r/a 
when  the  plasma  annulus  becomes  a  cylinder.  There  is  a  singularity  here,  however,  in  that 
the  speed  of  the  inner  radius  6  =  v(b)  =  u  does  not  smoothly  go  into  v(0)  =  0  as  6  —  0;  this 
will  require  special  treatment  later  in  the  discussion  of  the  implementation  of  these  equations 
in  the  code  and  the  concept  of  the  “therm  aiization*  procedure  (see  Section  III).  The  velocity 


*•  v  % 


wr! 
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profiles  v(r)  before  and  after  6  =  0  are  schematically  depicted  in  Fig.  2.  Assumptions  A1 
and  A2  allow  the  continuity  equation  (1)  to  determine  the  velocity  profile  (7)  in  terms  of 
the  outer  radius  a  and  its  velocity  u  (the  inner  radius  6  is  given  by  6(t)  =  a(t)  —  w,  where 
the  width  w  of  the  annulus  is  constant:  w  —  a  -  b  =  0). 


FIG.  2.  Velocity-  profiles,  (a)  v(r)  at  three  different  times  during  the  implosion  phase  as 
6  —  0,  showing  the  steepening  of  the  profile  near  r  =  b(t ).  (b)  v(r)  after  6  =  0  (pinch  and 
post-pinch  expansion  phase). 


Next,  the  force  equation  is  used  to  derive  an  equation  for  the  outer  radius  velocity  u, 


P  (j£  +  £ '  =  ~?p  +  |ixB. 


(8) 


Since  by  assumption  the  current  density  is  uniform  over  the  cross  section  of  the  skin,  the 

magnetic  held  B  =  B(r)6  is  easily  computed  to  be 

0  0  <  r  <  x  (inner  vacuum  and  core) 


B(r)  =  { 


2 1  r2  -  i2  . 

- 5 - r  x  <r  <a  (skin) 

cr  a3  —  x2 

21 

a  <  r  <  ry,  (vacuum) 


(9) 


v  cr 


Now  p,  t'(r),  J,  and  B(r)  are  all  known  quantities,  so  that  the  only  unknown  quantity  in 
(8)  is  the  pressure  profile  p(r).  Substitution  of  (2,7,9)  into  (8)  would  give  a  differential 
equation  (in  r)  for  the  pressure  in  terms  of  a(t),  u(t)  and  u(f).  However,  the  approach  used 
by  Guillory  and  Terry4  is  to  integrate  (8)  over  r  so  that  the  r-dependence  of  the  equation 
is  eliminated  along  with  the  Vp  term  (which  would  be  evaluated  at  the  endpoints,  r  —  b 
and  r  =  a).  This  then  gives  an  ordinary  differential  equation  in  time  for  the  outer  radius 
velocity  u  =  /(a.  b,  x.  u.  T).  This  equation  was  used  in  an  earlier  version  of  the  code  (prior 
to  March.  1985),  and  the  problem  with  it  seems  to  be  that  that  procedure  of  integrating 
over  r  gives  an  equation  which  is  inconsistent  with  energy  conservation.  This  will  not  be 
discussed  further  here,  but  it  will  be  shown  (see  Appendix  A)  that  the  present  method  does 
conserve  energy  (at  least  theoretically,  and  the  code  runs  seem  to  bear  this  out). 

In  order  to  eliminate  the  r-dependence  in  (8)  and  derive  an  equation  for  u(t).  first 
convert  (8)  into  an  equation  for  the  kintic  energy  density.  Multiplying  by  r  one  obtains 

p  ( Wt  +  V )  * t,S^r'  =  VP  ~  U  x  Sj-  (10) 

Multiplying  (1)  by  *r2  and  adding  the  respective  sides  of  that  equation  to  ( 10),  the  evolution 
of  the  kinetic  energy  density  k  —  kpv3  is  found  to  satisfy 


^^V.(Kt1)  =  -r.Tp-i£-(J7x;g). 
at  c 


HD 


Now  integrating  over  the  entire  volume  of  the  plasma,  the  total  kinetic  energy  A’(f)  obeys 
k(t)  =  -  J  <flr  v  ■  Vp  -  i  f  dV\  V  ■  U  *  B).  (12) 


£ 


The  integral  over  the  magnetic  force  term  is  only  over  the  volume  V,  of  the  skin.  In  obtaining 
the  left-hand  side  of  (12),  the  rule  for  time-varying  volumes  is  used 


4  f  dV*=  [  dV^-  -r  l  ds-vgK, 
Jvit)  Jvit)  at  Jsit) 


where  vs  is  the  local  velocity  of  the  surface  bounding  the  volume.  When  ( 13)  is  used  in  the 
volume  integral  of  (11),  the  term  arising  from  /  dVV  •  (kv)  =  f  ds  •  vk  cancels  the  similar 
surface  integral  term  in  (13)  if  the  local  surface  velocity  v$  is  equal  to  the  fluid  velocity  v 
at  that  point  on  the  surface.  Here,  of  course,  that  is  the  case,  and  the  conservative  form  of 
the  left-hand  side  of  (11,12)  merely  states  that  the  kinetic  energy  of  the  piasma  does  not 
change  due  to  particle  flux  through  the  bounding  surfaces  (at  r  =  6  and  r  =  a). 

The  pressure  term  in  (12)  can  be  converted  to  the  following  form 

-  J  dVv-  Vp  =  -  JdV’V-(pr)*  J  dV'p(V-e) 

=  —  2  ds  •  pv  (V  •  u)  I  dV  v 

Is  ~  J  (14) 

*  (V  •  u)  [J  <&c  P  -  J  pj 

=  (V  •  v)(Vepe  +  v,p,). 

In  the  second  line,  the  fact  that  in  this  model  the  divergence  of  the  velocity  field  is  indepen¬ 
dent  of  position  (6)  allows  it  to  be  moved  outside  the  integral.  In  the  third  line,  the  surface 
contribution  of  pv  vanishes  because  the  pressure  is  assumed  to  vanish  at  r  =  a  and  either 
the  pressure  is  zero  at  r  =  6  (when  b  >  0)  or,  when  6  =  0,  the  only  bounding  surface  is  at 
r  =  a.  Also  in  the  third  line,  the  integral  over  the  volume  has  been  broken  into  the  separate 
contnoutions  from  the  core  and  the  skin.  Finally,  in  the  fourth  line,  these  separate  contri¬ 
butions  are  expressed  in  terms  of  the  average  pressures  of  both  regions.  pe  S  V~l  J  dl'ep. 
etc.  The  separate  average  temperatures  and  total  number  of  particles  (ions  and  electrons) 


in  both  regions  are 


—  !'e,*tt(^e,i  ”  1  }Ic,»  —  c.tlc.ai 


where  a  third  assumption  is  made: 

A3.  The  ions  and  electrons  in  each  region  have  the  same  temperature:  T£  a  =  Tf  t  =  Te,s. 


Combining  (12,14,15),  the  equation  for  the  total  kinetic  energy  is 


k(t)  =  -—-(A refe  +  N.f.)  +  (vJB). 

a  •+■  o 


Here,  the  formula  for  V  -v  from  (6)  has  been  used,  and  the  magnetic  force  term  from  (12)  has 
been  abbreviated  by  ( vJB).  This  is  also  the  desired  equation  of  motion  for  the  outer  radius 
velocity  u(t):  substituting  v(r)  from  (7)  into  the  definition  of  the  kinetic  energy  density, 
integrating  over  the  total  volume,  then  differentiating  with  respect  to  time  and  using  this 
for  the  left-hand  side  of  ( 16 )  produces  an  equation  of  the  form 

Ciii  —  —  C?u'  +  C^cTc  +  C$tTa  —  C4/2.  (1<) 

Since  J  and  H(r)  are  also  known,  the  ( vJB )  term  can  be  evaluated  to  give  the  J2  term  in 
(17).  The  coefficients  Cx  are  functions  of  a,  b.  z.  Zc,  Z,  and  are  explicitly  given  in  Appendix 
B.  This  equation  is  coupled  to  the  following  temperature  equations  and  should  also  be 
accompanied  by  a  fourth  first-order  equation,  a(t)  =  u{t). 

The  derivation  of  the  equations  governing  the  thermal  energy  or  temperature  in  each 
region  (skin  and  core)  begins  with  the  pressure  equation  found  in  Braginskii5 

|  [^  -5-  V  -  (ps)]  =-pV-v-r,J2-R.  (18) 

Here,  the  resistivity  rj  is  e~x  and  R  is  the  radiated  power  density.  From  (8)  recall  that 
because  the  velocity  and  magnetic  fields  depend  on  r,  so  must  the  pressure  (and  in  fact. 
(8)  should  probably  be  interpreted  as  an  equation  for  the  pressure  profile).  However,  as 
in  the  case  of  the  kinetic  energy  equation,  we  will  integrate  (IS)  over  the  volume  in  order 
to  find  an  equation  for  the  average  pressures  (or  temperatures)  of  each  region.  Therefore, 
integrating  (IS)  over  the  core  and  skin  separately,  one  obtains 

|  ;/ ^  -T  /ffv;  V  .  (pv)  =-J  dV\p(T  .  V)  -  j  dV'R ,  (19a) 

l !  /  ^  ~  J  v-(p^  =  -  f  p(v  •  £)  -  /  ^  R-  (19b) 

Now  using  the  same  rule  for  time- varying  volumes  as  expressed  in  (13).  these  equations 
become 


T 


(20b) 


I '  K"  MW<1W  i«  hiji  .ujrm^hmmuui 

\~  f  &'.?  =  - 1  dv.P{v-v)  + 1  dv.  t,j3  -  j  dv.R 

This  step  requires  one  important  subtlety,  however.  In  using  the  rule  (13),  one  expects  the 
surface  integral  term  from  (13)  involving  the  local  boundary  velocity  to  cancel  the  surface 
integral  arising  from  the  divergence  theorem  applied  to  the  V  •  (pv)  term  on  the  left-hand 
side  of  (19),  which  involves  the  local  fluid  velocity.  Note  that  one  of  the  surfaces  involved 
in  the  surface  integrals  around  both  the  skin  and  the  core  is  the  surface  at  r  =  x  =  a  -  6, 
the  boundary  between  the  skin  and  the  core.  For  these  surface  integrals  to  cancel  then, 
one  must  require  that  this  boundary  move  with  the  local  fluid  velocity:  i  =  r(x).  where 
v(x)  is  the  velocity  profile  (7)  evaluated  at  r  =  x.  This  in  turn  implies  that  the  skin  depth 
must  vary  a s6  —  a  —  x  =  u-  u(x),  as  opposed  to  obeying  any  other  phenomenological  rule 
(such  as  6  —  ( c/ue ).  or  6  =  constant).  This  is  just  a  statement  that  thermal  energy  does 
not  diffuse  from  one  region  to  the  other,  as  was  originally  assumed.  Indeed,  integrating  the 
continuity  equation  (1)  over  the  volumes  of  the  skin  and  core  separately,  one  finds  that  such 
a  requirement  on  the  motion  of  the  boundary  between  them  simply  means  that  the  total 
number  of  particles  in  each  region  remains  constant:*  thus,  thermal  energy  is  not  diffused 
by  particles  of  the  higher  temperature  in  the  skin  moving  into  the  lower  temperature  core 
region.  If  such  diffusion  were  to  be  allowed,  one  would  have  to  include  other  terms  in  the 
pressure  equation  (18).  As  it  stands,  however,  this  requirement  on  the  variation  of  the  skin 
depth  is  necessary  for  total  energy  conservation. 

With  (15).  these  equations  can  be  converted  into  equations  for  the  average  temperatures 
of  each  region.  In  order  to  include  the  energy  required  for  ionization,  however,  it  is  useful 
to  introduce  an  alternate  temperature  which  measures  the  total  internal  energy  Q  of  the 
electron  fluid.  As  defined  by  Terry  and  Guillory,6  this  is 

<«> 

*  Actually,  this  integration  would  show  that  the  total  number  of  ions  in  each  region  is 
constant:  due  to  the  provision  in  the  model  for  ionization  depending  on  temperature,  the 
number  of  electrons  in  each  region  may  vary. 


Vi.' 
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where  &{Q)  is  the  specific  chemical  potential,  the  energy  required  to  ionize  an  atom  to  a 
charge  Z  while  maintaining  an  ambient  electron  temperature  Tt.  When  the  quantities  Q 
and  Z  are  calculated  properly  in  a  collisional  radiative  equilibrium  model,  there  will  always 
be  a  one-to-one  correspondence  between  Q  and  Te.  The  inverse  functions  T{Q)  and  Z{T) 
then  provide  an  equation  of  state,  and  only  the  time  development  of  Q  need  be  computed. 
Using  (15),  the  spatially  averaged  pressures  in  (20a.b)  can  be  recast  in  terms  of  the  total 
internal  energy  of  the  fluid 


2^'c,»Pc,t  —  Uc, j  —  c,»!A,jQc,j  T  Te><j. 


(22) 


Since  t  =  nV'CiJ  is  time  invariant  and  (V-r)  is  uniform  in  this  model,  the  time  development 
of  Qca  is  given  by 

l\e±(ZcQc  -r  fc)  =  ~Ne(Zc  +  1  )fc«?e)(V  •  v)  -  Vc5c,  (23a) 

§ Ns^(Z.Qt  -r  t.)  =  -N,{Z,  +  l)f.(&,)(V  •  v)  -  V0r)J‘  -  V,R..  (23b) 

Here  we  have  used  pV  =  .V* ( Z  +  1)T  on  the  right-hand  sides  (instead  of  (22))  because  only 
the  internal  energy  of  the  free  electron  fluid  contributes  to  the  p  dV  work.  This  modification 
of  Braginskii’s  heating  relation  involves  only  the  partition  of  thermal  energy:  the  source 
terms  are  undisturbed.  Thus,  the  volume  averaged  radiated  power  density  Re,,  has  been 
defined,  and  the  fact  that  the  current  density  J  is  uniform  over  the  skin  in  this  model  has 
also  been  used. 

The  time  derivative  on  the  left-hand  side  of  (23)  becomes 

U  =  \N\'ZQ^  Z£)  +  T), 

(24) 

=  \S\2'Q^2  +  T')Q, 


where  Z'  and  T'  are  derivatives  with  respect  to  Q.  In  principle,  this  form  for  V  could  be 
numerically  implemented  in  the  code,  but  instead  the  following  expressions  for  U  and  V 
have  been  used: 


uc.s  =  §-V’.,(Zc.,  +  i)gc<# 

=  l*TlcJZc,.  +  i)Qe,a 


(25) 


This  assumes 


A4.  The  energy  of  an  ion  is  §<3  instead  of  §T.  This  is  merely  a  deSnition:  for  purposes  of 
total  energy  inventory,  as  long  as  this  deSnition  is  consistently  used  then  it  should  not 
affect  the  conservation  of  energy. 

The  expression  for  U  in  (25)  is  a  good  approximation  to  the  exact  form  (24)  in  the  high 
temperature  limit  (where  Z'  %  0  and  f '  «  1).  For  this  approximation  we  therefore  assume 
A5.  The  rate  of  change  of  the  charge  state  over  a  single  timestep  is  negligible 
If  assumption  A5  is  violated  it  could  make  a  difference  in  determining  if  energy  is  conserved. 
If  the  charge  state  does  change  rapidly,  one  should  use  an  iterative  method  to  conserve  energy 
at  each  timestep  (such  as  that  used  in  the  thermalization  step  described  in  Section  IV). 

With  (25)  substituted  for  the  left-hand  sides  of  (23),  the  energy  equations  in  the  core 
and  skin  become 

*•— sirs <*■> 

<5.=  (=6b) 

o  a  v  o 

Here,  we  have  inserted  the  definition  of  (V  •  v)  from  (6),  and  the  coefficients  a  are  aCi,  s 
[§n(Ze.,  rl)]"1.  These  equations  close  the  set  of  ordinary  differential  equations.  As  stated 
earlier,  they  are  to  be  solved  in  conjuction  with  rules  for  obtaining  T(Q),  Z(T)  and  the 
radiated  power  density  R(T,  n,  Z). 


Ill  THERMALIZATION  of  KINETIC  ENERGY 


This  Section  discusses  the  procedure  which  is  used  in  the  code  to  convert  some  of  the 
kinetic  energy  of  the  plasma  at  assembly  time  (when  6  — »  0)  into  thermal  energy.  There  are 
two  reasons  for  making  such  a  conversion: 

1)  In  the  real  physical  system,  the  kinetic  energy  gained  during  the  implosion  (or  run-in) 
phase  is  partially  converted  into  thermal  energy  as  the  various  innermost  components 
of  the  gas-puff  ‘‘collide”  on  axis.  The  actual  physical  mechanism  for  this  rapid  heat¬ 
ing  probably  involves  collisional  effects,  counterstreaming  plasma  instabilities,  possibly 
shock  heating  and  other  complicated  processes.  The  theory  contained  in  this  report 
and  implemented  in  the  code,  however,  cannot  (and  does  not  intend  to)  treat  these 
physical  effects  in  any  real  way. 

2)  In  the  equations  of  motion  derived  in  Section  II  and  explicitly  displayed  in  Appendix 
B.  the  velocity  profile  v(r)  changes  abruptly  when  6  =  0  (see  Fig.  2):  before  6  =  0,  fluid 
elements  at  the  inner  radius  6  move  inward  with  speed  u  =  r(6)  =  v{a)  (see  Section 
II,  assumption  A.2)  which,  of  course,  becomes  very  large.  After  6  =  0,  however,  the 
fluid  velocity  at  r  =  6  =  0  is  suddenly  zero:  the  fluid  elements  approaching  the  axis  are 
(artificially)  abruptly  stopped.  Figure  3  shows  a  schematic  picture  of  the  kinetic  energy 
profiles  (~  r2(r).  from  Fig.  2)  before  and  after  6  =  0.  Evidently,  the  total  kinetic  energy 
(the  integral  under  these  curves)  contained  in  the  plasma  is  greater  when  6  ;>  0  than 
when  6  =  0.  due  to  the  abrupt  change  in  velocity  of  the  fluid  elements  near  r  =  6  %  0. 
This  discrepency  in  kinetic  energy  must  be  corrected  for,  or  else  the  plasma  w'ill  continue 
to  evolve  (with  the  new  force  equation  coefficients  valid  for  6  =  0,  Eqs.(B.lO))  with  less 
total  energy  than  it  had  before  the  assembly  time:  this  suggests  the  natural  solution 
described  below  for  converting  this  excess  kinetic  energy  of  fluid  elements  near  the  axis 
into  thermal  energy  of  the  entire  core. 


The  conversion,  or  thermalization  process  used  in  the  code  proceeds  as  follows:  assume 
that  at  the  beginning  of  timestep  n„  (time  ta  =  naAt),  the  inner  radius  6  is  for  the  first 
time  less  than  some  arbitrary  small  critical  radius  r0  (usually  chosen  to  be  1mm).  This  is 
called  the  assembly  time,  or  thermalization  timestep,  and  this  will  be  the  step  in  which  the 
equations  of  motion  will  be  changed  from  those  that  hold  before  6  =  0  (Eqs.(B.14)  with 
coefficients  (B.9))  to  those  valid  after  6  =  0  (Eqs.(B.14)  with  (B.10)).  Furthermore,  after 
the  calculations  using  the  present  value  of  6  are  finished,  6  will  be  set  to  zero. 


The  kinetic  energy  in  the  velocity  profiles  in  Figs.(2,3)  is  computed  in  (B.3).  The 
difference  in  this  kinetic  energy  (to  be  thermalized)  is  simply 

2o6  4a262  ln(a/6)  ' 


IK  =  \nlu7 


(27) 


_(a-i-6)2  (a  -  6)(a  6)3  J  ' 

This  energy  is  to  be  distributed  over  all  the  particles  (ions  and  electrons)  in  the  core  as  an 
increase  in  thermal  energy  A Q: 


A  K  =  §.VeAQ  =  §(Ze  +  1)JV‘A  Q. 


(28) 


Here.  .V*  =  (^£/m,)[(x2  -  b~)/(a2  -  62)!  is  the  number  of  ions  in  the  core  (a  constant,  see 
Section  II).  and  Ne  =  (Ze  —  1)A*C‘  is  the  total  number  of  electrons  and  ions  in  the  core  at 
the  beginning  of  this  timestep.  From  (28),  the  new  thermal  energy  per  particle  in  the  core 
after  thermalization  would  be 


Q[l)  =  A<?  +  Qi0)  = 


2A  K 


3  (Zj0)  +  1)JV* 


-Q 


(0) 


(29) 


Now.  Q(c0)  and  Zi0)  are  the  values  of  the  thermal  energy  per  core  particle  and  average  core 
ion  charge  state  at  the  beginning  of  the  timestep.  The  problem  with  using  Qi1''  for  the 
new  thermal  energy  in  the  rest  of  the  timestep  is  that  it  is  generally  too  large  and  will  not 
conserve  energy.  To  see  this,  consider  the  total  amount  of  energy  £0  available  for  thermal 
energy  at  the  beginning  of  the  timestep: 


£0  =  A  A'  -  §(Zj0)  +  1)A T*qW.  (30) 

Using  the  equation  of  state  (see  Section  IV),  the  new  thermal  energy  per  core  particle  QcX  i 
would  produce  a  core  temperature  Tlli(Q{‘'1)  and  new  ion  charge  state  Zi*'(7'e1))  which 
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would  be  much  greater  than  their  values  T<  ,  Zr‘  at  the  beginning  of  the  timestep.  Thus, 
the  new  total  core  thermal  energy  would  be 


LtW  =  §.y"Z<»  + 


(31) 


z\ 


U) 


Zi0}  + 1 


E0. 


For  neon,  the  typical  case  is  that  zi°*  %  1,  while  the  initial  thermalization  estimate  is 


(i) 


enough  to  strip  all  but  the  innermost  shell  of  electrons  giving  zi1^  %  8.  Therefore,  if  Qc 


>(»> 


were  used,  the  total  core  thermal  energy  Ui1'*  would  be  about  4-5  times  the  total  amount 
of  available  energy  E0  at  the  beginning  of  the  timestep.  Obviously,  what  has  happened  is 
that  this  energy  Eo  was  not  appropriately  split  into  thermal  and  ionization  energies. 

The  correct  way  to  implement  the  thermalization  would  be  to  solve  the  implicit  equation 
for  Q(c1] 


§^e‘[z<1>(Q(l))  +  l]  Qil)=E0,  (32) 

using  the  equation  of  state  (see  Section  IV).  That  is,  one  must  adjust  the  new  value  of  Qc  so 
that  the  corresponding  number  of  electrons  Zc  per  ion  have  the  self-consistent  thermal  energy 
ZeQe  to  give  a  total  core  thermal  energy  Uc  equal  to  the  available  energy  Eo-  Therefore,  at 
the  beginning  of  the  thermalization  timestep,  the  main  routine  IMPLODE  calls  the  subroutine 
THERM,  which  uses  an  iterative  procedure  to  determine  the  new  Qe.  This  is  described  briefly 
below  and  illustrated  in  Fig.  4.  Since  is  too  large,  a  new  trial  value  is  taken  to  be 
the  average  of  the  initial  value  Q^°'  and  Qi^.  The  total  core  thermal  energy  is  then 
computed  with  the  new  charge  state  Z'c2]  and  compared  with  the  consant  value  Eo-  If  Crj2) 
is  still  greater  than  Eo,  Qi is  chosen  as  an  averge  between  Qe°'  and  if  Cre2)  <  Eo, 
Q(c3)  is  taken  to  be  the  average  of  Qi1^  (too  high)  and  This  method  of  “splitting  the 
difference”  continues  until  the  values  of  and  Eo  agree  to  within  some  tolerance  (usually 
picked  to  be  about  10J). 


E0  =  AA'  +  §(Ze(0)  +  1)-V*Q<°) 


Qix)  =  A<5  +  <?<0)  = 


2AA 


3(zi0)  +  i)i\r* 


+  Qi0) 


<?m.»  =  Ql°] 
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FIG.  4.  Flowchart  illustrating  therm  alizatioD  procedure  used  in  subroutine  THERM 


As  an  example  of  the  numbers  involved,  in  a  particular  code  run  the  core  plasma  had 
the  following  parameters  at  the  beginning  of  the  thermalization  timestep:  A'*  ss  3.5  x  1018, 
=  20.1  eV,  Te(0)  =  4.5  eV,  Zi°'  =  1.38.  With  a  kinetic  energy  A K  =  1795  J  to  be 
thermalized,  these  produce  the  initial  values  A Q  —  884  eV,  ui°^  —  40.1  J,  and  Eo  —  1835  J. 
After  n  =  10  iterations,  the  final  thermalized  values  were  determined  to  be  Qin>  =  273  eV, 
Tin)  =  158  eV,  zjn)  =  8.0,  and  LTin)  =  1833  J. 

IV  ATOMIC  PHYSICS:  IONIZATION  AND  RADIATION 

The  reason  for  developing  this  model  of  an  imploding  plasma  radiation  source  is  to 
be  able  to  compute  the  amount  of  radiation  produced  by  the  plasma.  The  description  of 
the  atomic  physics  involved  in  collisions!  ionization  and  radiation  can  be  on  many  different 
levels  and  can  be  based  on  several  different  models.  The  primary  concern,  however,  is  to 
keep  this  simulation  as  sip’ple  as  possible  so  that  in  can  be  used  in  conjunction  with  a 
simulation  of  the  entire  puised-power  device  for  purposes  of  design  and  operating  parameter 
optimization.  Thus,  for  purposes  of  determining  ionization  levels  and  radiation  we  have 
used  the  results  of  Terry  and  Guillory:6  these  are  in  the  form  of  piecewise  numerically-fit 
curves  giving  an  equation  of  state,  and  average  ion  charge  state  and  radiated  power  density 
as  functions  of  plasma  electron  temperature  and  density  for  neon  and  argon.  The  original 
data  for  the  equation  of  state  and  radiated  power  density  were  obtained  from  collisional 
radiative  equilibrium  calculations  due  to  D.  Duston. 

A  EQUATION  OF  STATE  and  AVERAGE  ION  CHARGE  STATE 

An  equation  of  state  is  used  to  determine  the  branching  of  ‘Thermal  energy”  per  electron 
§<?  into  actual  thermal  energy  §  T  and  the  energy  ©  needed  to  ionize  the  average  atom  to  a 
change  state  consistent  with  the  electron  temperature  T.  As  in  equation  (21),  this  quantity 
Q  has  been  defined  as 

\Z{T)Q  =  \Z{T)T{Q)  +  ©(D,  (33) 
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where  0(T(Q))  is  the  energy  required  to  ionize  an  atom  to  an  average  charge  state  of 
Z{T(Q)).  When  Q  is  incremented  at  each  timestep,  the  subroutine  TFUNC  is  called  to 
determine  T(Q)  and  subroutine  ZFUNC  is  then  called  to  determine  Z(T).  These  values  of 
T  and  Z  are  used  in  computing  the  coefficients  in  the  differential  equations  for  the  next 
timestep.  For  neon,  a  numerical  fit6  for  the  branching  ratio  B(Q)  is  used  where 

B(Q)  =  -9—  -  1  _ — ~21Q1 — .  (34) 

yW'  T{Q)  3Z(Q)T{Q)  1  1 

Thus,  knowing  Q.  TFUNC  simply  computes  B(Q)  and  hence  T{Q)  =  Q/B{Q).  (Actually, 
the  piecewise  fit  in  TFU.VC  is  based  on  the  logarithms  of  Q  and  B. )  The  graph  of  In  T  as  a 
function  of  In  Q  for  neon  is  shown  in  Fig.  5.  Subroutine  ZFUNC  uses  a  numerical  piecewise- 
iinear  fit6  to  a  curve  Z(T ),  as  shown  in  Fig.  6.  The  code  is  also  capable  of  treating  an  argon 
plasma,  but  results  will  only  be  presented  for  neon. 


B  RADIATED  POWER  DENSITY 

la  order  to  compute  the  radiated  power  density  from  the  imploding  plasma,  we  have 
used  the  results  of  more  sophisticated  atomic  physics  computer  calculations  produced  by 
Terry.6  These  are  in  the  form  of  piecewise-spline  fitted  numerical  tables  for  continuum  and 
line  radiation  as  a  function  of  electron  temperature  and  ion  density.  Four  types  of  radiation 
are  considered:  Bremsstrahiung,  free-bound,  and  A'-sheil  and  A-shell  X-rays.  The  radiated 
power  density  (W/cm3)  curves  for  these  types  of  radiation  for  neon  at  an  ion  density  of 
r.,  =  10 18  are  shown  in  Fig.  7.  Of  particular  experimental  interest  are  the  A-shell  emission 
in  neon  and  the  L-snell  line  in  argon:  these  Figures  show  that  both  of  these  lines  peak  at 
about  400  eV  for  the  respective  elements.  Each  type  of  radiation  is  produced  with  a  power 
density  that  scales  roughly  as  the  square  of  the  ion  density.  These  curves  are  the  result  of 
the  subroutine  CREMIT  obtained  from  R.  E.  Terry.6 

The  comparison  of  numerically  computed  radiation  yields  with  those  observed  in  exper¬ 
iments  are  based  on  several  aspects  of  the  radiation  pulse:  total  amount  of  radiated  energy 
f y V .  including  all  types  of  radiation),  total  radiated  X-shell  and  I-shell  X-ray  energy  (Yx 
and  Yi).  X-ray  efficiency  (ratio  of  YK  or  Yi  to  Vj)  ana  the  length  of  the  radiation  pulse.  We 
have  found  that  on  a  typical  code  run  with  plasma  parameters  which  can  be  matched  with 
a  corresponding  experiment,  the  radiation  yields  Y?  and  Yx  or  V'x,  do  not  agree  with  the 
experimentally  observed  values.  However,  as  the  plasma  parameters  are  varied  in  both  the 
numerical  simulations  and  the  experiment,  the  behavior  of  the  numerically  computed  yields 
appears  to  follow  that  of  the  experiment.  For  example,  fixing  the  initial  plasma  radius  (or 
nozzle  radius)  and  generator  voltage,  the  radiated  yields  increase  to  a  maximum  and  then 
decrease  again  as  the  mass  per  unit  length  m  is  varied:  this  is  observed  in  both  simulation 
and  experiment,  although  the  magnitudes  of  the  yields  do  not  agree  (the  simulation  values 
are  generally  greater). 


Because  of  this  similar  behavior,  we  have  attempted  to  bring  the  computed  yields 
into  better  agreement  with  the  experimentally  measured  values  by  introducing  an  overall 
radiation  multiplicative  factor  called  the  opacity .  The  way  in  which  this  artificial  factor 
is  set  (usually  a  number  between  zero  and  one)  is  described  in  the  next  Section;  once 
chosen  however,  (at  a  particular  set  of  parameter  values)  it  is  found  that  the  simulations 
are  in  much  better  agreement  with  experiment  over  a  wide  range  of  parameter  values.  The 
implementation  of  this  factor  is  simple:  at  each  timestep.  the  radiated  power  densities  Rc 
and  R,  must  be  computed  not  only  for  contributing  to  the  total  radiated  energy,  but  also 
for  acting  as  an  energy  sink  in  the  evolution  equations  (26)  for  Qc  and  Qa.  Thus,  the  main 
routine  IMPLODE  calls  the  subroutine  PLOSS.  which  in  turn  calls  CREMIT.  When  CREMIT 
returns  the  values  of  the  four  types  of  radiated  power  densities  for  the  current  temperature 
and  plasma  density.  P LOSS  simply  multiplies  each  by  the  value  of  the  opacity  factor  and 
and  returns  the  results  to  IMPLODE  for  use  in  the  evolution  equations  and  radiated  energy 
compilation.  Note  that  while  the  use  of  the  opacity  factor  usually  means  that  less  total 
energy  is  radiated,  it  also  means  that  the  radiation  is  less  of  a  sink  for  the  thermal  energy 
Q.  This  tends  to  keep  the  temperature  of  the  plasma  up.  but  actually  this  is  a  complicated 
effect:  not  only  is  the  radiation  a  temperature  dependent  quantity,  but  as  the  temperature 
remains  high,  so  does  the  plasma  pressure  during  the  pinch  phase  so  that  the  minimum 
radius  acheived  is  greater  (lower  pinch  density),  again  reducing  the  radiation. 


As  will  be  seen  in  the  next  Section  on  numerical  results,  the  peak  core  temperature 
typically  observed  at  pinch  time  (time  of  minimum  radius  or  maximum  density)  is  around 
400-500  eV.  compared  with  experimentally  measured  values  of  about  100-200  eV.  The  total 
radiations  yields,  however,  agree  well  with  the  experiment  (with  the  use  of  the  opacity 
factor).  This  is  somewhat  of  a  contradiction,  because  for  maximum  A'-shell  yield  (in  neon), 
the  radiation  curves  in  Fig.  7  indicate  that  at  pinch  the  temperature  should  be  precisely 
in  this  range,  yet  the  experiments  suggest  that  the  optimum  temperature  for  Yk  should 
be  nearer  100  eV.  Thus,  not  only  should  the  temperatures  in  the  simulation  remain  lower. 
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but  there  should  be  a  corresponding  shift  in  the  radiation  curves  to  lower  temperatures 
(otherwise,  very  little  radiation  would  be  computed).  One  possibility  is  that  the  equations 
of  state  that  are  used  are  not  appropriate:  that  is,  it  may  take  much  more  energy  than 
©  =  \Z{Q  —  T )  given  by  our  equations  of  state  to  ionize  an  atom  to  the  average  value  Z. 
Correcting  for  this  would  mean  dumping  more  available  “thermal  energy*  Q  into  ionization 
energy  than  presently  used:  this  would  act  to  keep  the  temperature  down.  The  effect  of  true 
plasma  opacity  to  X-shell  radiation  at  lower  temperatures  (in  neon)  on  the  ff-shell  yield 
has  also  been  suggested  ‘  as  a  phenomenon  not  accounted  for  in  the  radiation  curves  which 
could  possibly  lower  the  optimum  A'-shell  radiation  temperature  (in  neon)  to  the  100-200 eV 
level. 


V  NUMERICAL  RESULTS 


This  Section  discusses  the  use  of  the  code  IMPLODE  for  comparison  with  experiment  and 


parameter  optimization,  and  then  we  discuss  some  of  the  numerical  results.  First,  the  input 

parameters  available  to  the 

user  which  set  the  properties  of  gas-puff  and  its  environment 

and  which  remain  constant  through  a  single  code  run 

are 

Parameter 

j  Computer  1 

1  Variable 

1 

Definition 

Units 

j  Typical  j 

!  Value  i 

a0 

aO  j 

initial  radius 

cm 

i  1.25  cm  | 

da  | 

mass/length 

g/cm 

3.0  x  10-5  g/cm  i 

l 

Xl  1 

length 

cm 

4.0  cm  ! 

w 

w 

annulus  width 

cm 

0.3  cm 

A 

aion 

atomic  weight  i 

1 

20.183 

opacity 

opac 

opacity  1 

1 

0.25 

i>o 

deltO  i 

initial  skin  depth  i 

cm 

i  0.1  cm 

»-o 

r0  | 

critical  radius 

cm 

0.1  cm 

rw 

nr  ! 

wall  radius  i 

cm 

3.5  cm 

VG 

▼peak  i 

generator  voltage  j 

MV 

3.0  MV 

vo  r0  ;  initial  velocity  I  cm/ns  0.0 


The  values  given  are  for  a  typical  code  run  with  neon  in  the  GAMBLE  U  configuration 
of  Fig.  8  near  optimum  mass/length  at  an  open  circuit  generator  voltage  of  3MV  (peak  load 
current  1.2MA).  All  of  these  input  parameters  are  local  to  IMPLODE  and  its  accompanying 
subroutines,  except  for  VG  which  is  actually  input  with  the  system  configuration  parameters 
(the  remainder  of  which,  such  as  line  inductances,  are  not  listed  here).  The  value  of  the 
atomic  weight  A  should  be  either  20.183  (neon)  or  40.0  (argon).  In  principle,  other  species 
can  be  used;  presently,  however,  models  exist  in  the  code  for  only  neon  and  argon.  The 
initial  skin  depth  <5o  is  of  course  not  an  input  parameter  in  the  version  of  the  code  which 
does  not  have  a  skin  (see  Appendix  B).  The  critical  radius  r0  is  that  which  the  inner  radius 
6  of  the  gas-puff  must  reach  for  the  thermalization  step  to  occur.  The  parameters  of  most 
interest  for  design  optimization  are  ao,  p  and  VG.  Those  which  should  be  varied  to  analyze 
the  behavior  of  the  model  itself  are  w,  opacity ,  6q  and  r0.  The  parameters  £,  rw  and  A 
depend  on  the  experimental  situation,  and  »o  is  generally  always  taken  to  be  zero. 

A  setting  parameters  to  compare  with  experiment 

In  order  to  predict  optimal  performance  from  numerical  simulations,  one  must  first 
check  whether  the  code  results  agree  with  present  experimental  observations,  and  make 
necessary  adjustments  if  they  do  not.  The  first  step  in  setting  the  code  parameters  to 
coincide  with  experimental  parameters  is  to  choose  the  correct  system  configuration  in 
the  transmission  iine  code8  (TLC),  including  the  value  of  the  peak  generator  open  circuit 
voltage  VG.  The  comparison  with  experiment  is  best  made  in  terms  of  the  profile  of  the 
current  flowing  in  the  load  (the  PRS):  the  primary  points  of  comparison  are  the  peak  load 
current  Ip  and  the  time  tp  at  which  the  peak  occurs  relative  to  the  generator  discharge. 
This  comparison  can  be  made  with  almost  any  reasonable  values  of  the  other  parameters 
(although  one  should  try  to  get  an  estimate  of  the  appropriate  values  from  the  experiment), 
as  Ip  is  fairiy  insensitive  to  the  other  parameters  (a  slight  variation  in  Ip  with  p  at  fixed 
V'G  has  been  observed). 
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FIG.  8.  Transmission  line  configuration  for  GAMBLE  II  experiment  used  in  simulation. 

The  parameters  ao  and  p  are  fairly  difficult  to  determine  from  a  particular  experimen¬ 
tal  run.  Streak  photographs  of  visible  light  and  X-rays  from  the  implosion  and  data  from 
the  nozzle  manufacturer  show  that  the  gas-puff  is  not  axially  uniform:  it  tends  to  expand 
radially  as  it  moves  away  from  the  nozzle,  so  that  ao  and  w  appear  to  increase  along  the 
axial  direction.  Furthermore,  the  initial  plasma  density  profile  is  not  uniform  in  the  radial 
direction  across  the  annular  region.  These  effects  make  it  difficult  to  define  what  the  appro¬ 
priate  'outer  radius”  and  thickness  w  should  be  for  the  model.  The  experimental  measure  of 
the  mass  per  unit  length  p  of  the  gas  is  in  terms  of  the  initial  plenum  pressure.  It  is  believed 
that  over  the  parameter  range  studied  this  plenum  pressure  should  be  linearly  related  to  p, 
but  there  is  not  an  d  priori  model  for  this  relationship.  Moreover,  one-dimensional  simula¬ 
tions  which  include  the  effects  of  anomalous  resistivity9  in  the  early  phase  of  the  implosion 
suggest  that  some  mass  is  left  behind  as  the  load  is  snowplowed  up. 

It  will  be  shown  below  that  at  fixed  ao  and  Vc  there  is  an  optimum  value  of  p  in 
the  simulation  and  an  optimum  plenum  pressure  in  the  experiment.  Having  matched  a0 
and  V'a  to  the  experiment,  one  can  identify  the  optimum  value  of  p  with  the  optimum 
plenum  pressure.  This  value  of  p  must  be  interpreted  carefully:  it  may  measure  only  the 
mass  remaining  in  front  of  the  magnetic  piston  after  the  early  resistive  field  penetration 
has  ceased.  Whatever  the  real  meaning,  keeping  the  nozzle  radius  ao  constant,  the  optimal 
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values  of  n  and  plenum  pressure  can  be  compared  for  different  values  of  the  peak  current 
Ip  (or  Vq)  to  determine  a  fairly  linear  relationship. 

B  SETTING  THE  OPACITY 

In  order  to  scale  the  radiation  yield  produced  in  the  simulation  (t.e.,  determine  the 
opacity),  the  following  procedure  is  used.  Matching  ao  and  Ip  (or  V'c)  to  experimental 
values,  the  opacity  is  set  to  unity  and  n  is  varied  to  compute  radiation  yields  Y7,  Yk  and  Yi 
as  a  function  of  p..  At  the  value  of  fi  at  which  Yk  is  maximum  for  neon  (Yi  for  argon),  we 
compare  Yk  (Yl)  with  the  experimental  value.  As  stated  in  Section  IV,  the  radiation  yields 
computed  in  the  simulation  with  unit  opacity  will  generally  be  higher  than  that  measured 
in  experiment.  Thus,  one  proceeds  to  reduce  the  opacity  and  compute  the  radiation  yields 
at  fixed  ao  and  V’c,  keeping  n  set  at  the  optimum  value.  Once  a  value  of  the  opacity  is 
found  that  produces  yields  that  are  in  good  agreement  with  experiment,  the  mass/length  is 
again  varied  (at  this  new  opacity)  in  the  neighborhood  of  the  optimum  n:  this  is  because  the 
optimum  a  is  generally  shifted  slightly  higher  as  the  opacity  is  reduced.  This  procedure  can 
be  repeated  until  the  desired  agreement  with  experiment  is  reached  for  some  combination 
of  opacity  at  optimum  (1.  For  neon,  with  nozzle  radius  ao  =  1.25  cm  and  Vc  =  3.0  MV 
(Ip  =  1.2  MA),  the  optimum  mass/length  was  found  to  be  ^  =  30/xg/cm  with  an  opacity 
of  0.25.  Ia  the  version  of  the  code  without  the  skin  (see  Appendix  B),  the  optimum  mass 
was  38^g/cm  with  an  opacity  of  0.12. 

C.  RADIATION  SCALING  WITH  CURRENT  AND  MASS/LENGTH 

It  is  hoped  that  setting  the  opacity  in  this  way  (so  that  the  computed  yields  agree 
with  experiment  at  a  single  data  point)  will  allow  one  to  vary  the  other  experimentally 
variable  parameters  (such  as  n  and  Vc)  and  obtain  similar  agreement.  If  this  is  indeed  the 
case,  then  one  would  feel  justified  in  varying  other  parameters  (such  as  ao)  which  are  not 
easily  changed  in  experiment  in  order  to  predict  better  experimental  design.  Therefore,  the 
following  series  of  simulations  were  performed  to  compare  with  similar  experimental  results. 


Fixing  ao  to  the  experimental  value  of  1.25  cm,  the  yield  was  investigated  at  eight  values 
of  Vc  (or  Ip)  with  neon  gas.  The  results  for  the  X-shell  X-ray  yield  Yk  with  respect  to 
variation  of  mass/length  at  each  of  the  eight  current  levels  is  shown  in  Fig.  9.  As  mentioned 
several  times  above,  there  is  an  optimum  value  of  p  for  each  current  level,  varying  between 
15  and  45  Mg/cm  (the  version  of  the  code  with  the  skin  was  used).  The  single  data  point 
which  was  used  to  normalize  the  opacity  to  match  experiment  is  the  optimal  point  on  the 
curve  with  Ip  =  1.23  MA.  Corresponding  experimental  results  for  five  current  levels  are 
shown  in  Fig.  10:  note  that,  again,  these  yields  are  experimentally  measured  with  respect 
to  the  initial  gas  nozzie  pressure.  The  general  trend  of  the  curves  here  compares  well  with 
the  simulation  results:  the  peak  yield  on  several  experimental  curves  agree  to  within  about 
109c  with  the  peaks  in  the  simulation  curves,  when  the  current  levels  can  be  matched.  This 
comparison  can  be  made  somewhat  clearer  by  plotting  the  simulation  results  of  Fig.  9  against 
pressure,  instead  of  If  one  plots  the  optimum  p.  against  optimum  pressure  for  five  of  the 
simulation  curves  (or  interpolations  of  them)  which  roughly  match  the  current  levels  of  the 
experimental  curves,  one  obtains  the  linear  relationship  shown  in  Fig.  11.  This  can  then  be 
used  to  translate  the  simulation  curves  into  functions  of  pressure,  as  shown  in  Fig.  12,  which 
now  can  be  directly  compared  with  the  experimental  results. 

Probably  the  most  notable  difference  is  the  relative  widths  of  the  curves:  the  exper¬ 
imentally  measured  radiated  energy  appears  to  be  slightly  less  sensitive  to  variations  in 
pressure  (or  n).  One  explanation  for  the  narrower  curves  produced  by  the  simulation  could 
be  that  the  radiated  energy  density  curves  used  (Fig.  7)  are  too  narrow  functions  of  temper¬ 
ature.  To  see  this,  one  first  needs  to  understand  why  there  is  an  optimum  n  for  fixed  values 
of  the  other  parameters.  A  simple  reason  based  on  the  thermalization  concept  described 
in  Section  III  is  as  follows:  If  one  assumes  that  the  evolution  of  the  plasma  annulus  during 
the  implosion  or  run-in  phase  (».e.,  before  assembly  time,  or  thermalization)  is  driven  by  a 
magnetic  force  which  is  constant  in  time,  then  the  kinetic  energy  imparted  to  the  plasma 
would  be  the  same  for  all  values  of  n-  This  is  because  the  constant  force  would  act  over 
the  same  radial  distance  (ao  -  (r0  -f  ic)),  even  though  the  time  required  for  run-in  would 
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FIG.  10.  Experimental  results  for  neon  A"-shell  radiated  energy  Yk  as  a  function  of 
plenum  pressure  for  five  values  of  peak  load  current  Ip  (labelling  the  curve  in  MA).  The 
nozzle  used  in  experiment  has  a  radius  of  1.25  cm. 
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FIG.  12.  Numerical  results  from  Fig.  9  plotted  in  terms  of  nozzle  pressure  using  the 
best-fit  linear  relationship  of  Fig.  11. 


This  is  where  the  mass/length  makes  a  difference:  the  larger  the  value  of  pi,  the  more  ions 
(and  electrons)  there  are  in  the  core,  so  the  thermalization  process  distributes  less  energy 
per  particle.  Thus,  for  small  values  of  pi.  the  core  temperature  increases  so  much  during 
thermalization  that  the  resuiting  pressure  bounces  the  plasma  back  out  before  it  reaches  a 
density  which  produces  a  significant  amount  of  radiation  (R  ~  n2).  Furthermore,  the  tem¬ 
perature  would  overshoot  the  optimum  value  for  the  X-ray  emission  On  the  other  hand,  if 
pi  is  very  large,  the  temperature  due  to  thermalization  is  too  small  for  the  plasma  to  ever 
reach  optimum  X-ray  emission  temperatures  at  pinch  time,  and  hence  only  a  small  amount 
of  radiation  is  observed.  Therefore,  somewhere  between  small  values  and  large  values  of  pi 
there  must  be  values  where  significantly  more  radiation  is  emitted. 

Both  the  real  physical  system  and  the  simple  model  described  in  this  paper  are  of 
course  more  complicated  than  the  explanation  given  above.  However,  what  is  true  is  that 
there  is  a  delicate  balance  between  the  mass/length,  the  temperature  achieved  at  pinch, 
the  density  achieved  at  pinch  and  the  radiation  produced:  the  interaction  of  these  effects 
are  responsible  for  the  existence  of  an  optimum  mass/length.  Because  of  this  temperature 
dependence  of  the  radiation  then,  the  off-optimum  performance  can  be  understood  also.  If 
the  radiation  curves  in  Fig.  7  were  extremely  narrow  in  temperature,  then  only  the  optimum 
pi  value  would  radiate  significantly  and  any  case  with  off-optimum  pi  would  produce  very 
little  radiation.  Thus,  the  widths  of  the  yield  curves  in  Figs.  9- 12  must  be  related  to  the 
widths  of  the  temperature  dependence  of  the  radiation  in  Fig.  7.  It  is  evident  then  that  the 
actual  physical  system  radiates  significantly  over  a  slightly  larger  temperature  range  than 
that  assumed  in  the  model. 

The  variation  of  A'-shell  radiated  energy  Yjc  with  peak  load  current  Ip  xa  Figs.  9- 12 
is  quite  dramatic.  In  Fig.  13.  we  plot  the  value  of  Yk  at  optimum  pi  for  each  current  level 
curve  in  Fig.  9.  This  log-log  plot  has  a  best  fit  slope  of  4.02,  which  implies  a  scaling  law 
Yk  ~  I4-  Furthermore,  plotting  in  Fig.  14  the  value  of  Yk  at  optimum  pi  against  pi  for  each 
curve  in  Fig.  9.  one  finds  a  Yk  ~  M3  scaling.  These  scaling  laws  agree  with  those  found  in 
experiment,  shown  in  Fig.  15. 
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FIG.  14.  For  each  constant-current  curve  in  Fig.  9,  the  value  of  Yk  at  optimum  n 
plotted  against  the  value  of  fi.  The  best-fit  slope  is  2.1. 


FIG .  15.  Experimental  scaling  laws  corresponding  to  those  in  Figs.13. 14  using  the  mea¬ 
sured  results  shown  in  Fig.  10. 
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Such  scaling  laws  have  been  explained  in  the  past  on  the  basis  of  very  simple  physical 
arguments.  These  arguments,  however,  do  not  take  into  account  the  temperature  depen¬ 
dence  of  the  radiation  and  thus  do  not  contain  the  concept  of  an  optimum  mass/length. 
In  a  greatly  simplified  version  of  the  theory  presented  in  Section  II,  the  radiation  can  be 
treated  as  a  perturbation  on  the  dynamics  of  the  implosion  after  the  assembly  (thermal- 
ization)  time,  and  used  to  calculate  curves  similar  to  those  in  Fig.  9.  In  that  theory,  the 
ideas  of  optimum  mass/ length  and  the  broadening  of  the  yield  curves  being  associated  with 
the  width  of  the  radiation  dependence  on  temperature  are  made  explicit.  In  addition  to 
producing  the  scaling  laws  observed  above,  this  theory  explains  the  presence  of  an  optimum 
mass  in  terms  of  the  idea  of  an  optimum  implosion  velocity  at  assembly  (thermalization) 
time:  that  is,  the  value  of  the  optimum  mass  for  a  given  current  is  that  which  reaches  a 
particular  velocity  (independent  of  current)  at  assembly  time.  Evidence  for  this  idea  has 
been  observed  in  experiment10  and  is  provided  by  the  numerical  results  shown  in  Fig.  16.  In 
(16a).  we  plot  the  velocity  at  thermaiization  ttr  of  the  optimum  mass/length  case  against 
the  load  current  from  Fig.  9;  in  (16b),  ut  for  the  optimum  mass  case  is  plotted  against  the 
value  of  fi.  That  the  value  of  uT  is  indeed  fairly  constant  in  these  optimum  cases,  however, 
may  only  indicate  that  the  simplified  treatment  which  predicted  such  a  property  succeeds 
in  replacing  the  more  complicated  equations  of  the  model  which  the  code  actually  evolves. 
The  existence  of  an  optimum  velocity  in  the  real  physical  system  is  still  an  open  question. 

D  EVOLUTION  OF  DYNAMICAL  VARIABLES 

As  an  example  of  the  numerical  solution  to  the  equations  of  motion  given  in  Section 
II  and  Appendix  B.  the  evolution  of  several  of  the  dynamical  variables  and  related  physical 
characteristics  of  the  imploding  plasma  for  a  particular  set  of  input  parameters  are  presented 
in  this  subsection.  The  parameters  are  in  fact  those  given  in  the  table  at  the  beginning  of 
this  Section,  except  for  which  was  taken  to  be  the  value  of  c/ae  for  the  initial  density 
(this  is  actually  much  too  small  for  a  skin  depth,  as  will  be  seen). 
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FIG.  16.  For  each  constant-current  curve  in  Fig.  9,  we  plot  the  velocity  u t  at  thermal 
ization  for  the  optimum  mass/length  case  as  a  function  of  (a)  the  value  of  Ip  on  the  curve 
and  (b)  the  value  of  n  for  the  optimum  case. 


Fig.  17(a)  shows  the  trajectory  of  the  outer  radius  a,  beginning  at  ao  =  1.25  cm,  ther- 
malizing  at  t  =  106ns  (when  b  =  r0  =  1mm,  a  =  4  mm),  pinching  at  t  =  117cm  and 
bouncing  back  out  (although  the  code  run  of  200  ns  is  not  long  enough  for  it  to  reach  the 
wall  radius.  rw  =  3.5  cm).  The  outer  radius  velocity  u(t)  is  shown  in  (b)  accelerating  to  a 
maximum  inward  speed  of  0.033  cm/ns  at  thermalization;  then,  with  thermalization  produc¬ 
ing  a  significant  plasma  temperature  and  pressure  (see  (d)),  the  inward  motion  immediately 
begins  to  slow,  passing  through  u  =  0  at  pinch  and  finally  reaching  a  fairly  constant  out¬ 
ward  speed  when  the  temperature  and  current  decrease  significantly.  The  density  rises  quite 
rapidly  to  a  peak  of  about  2.6  x  1019cm“3  at  pinch  in  (c).  The  core  temperature  history 
in  (d)  reveals  the  abrupt  effect  of  the  thermalization  procedure:  at  t  —  106  ns,  Te  rises 
from  less  than  5eV  to  158  eV  (the  details  of  the  thermalization  are  in  the  example  numbers 
used  at  the  end  of  Section  III).  With  this  boost.  Tc  continues  to  rise  due  to  compressional 
heating,  reaching  a  maximum  of  about  460  eV  near  pinch.  As  stated  earlier,  this  is  the 
optimum  n  case  for  the  present  current  level,  and  this  peak  temperature  is  in  the  correct 
range  for  maximum  A'-shell  X-ray  emission:  the  peak  is  reached  just  slightly  before  pinch, 
however,  because  as  the  temperature  and  density  become  optimum  for  radiation,  the  surge 
of  radiation  acts  to  cool  off  the  plasma.  This  radiation  pulse  is  shown  in  (e).  where  the  solid 
line  is  the  total  radiated  power,  the  dotted  line  (nearly  invisible  under  the  solid  line)  is  the 
A'-shell  A'-rav  radiated  power,  and  the  curve  marked  with  an  “L*  is  the  L-snell  radiated 
power.  Evidently,  since  the  pinch  temperature  is  appropriate  for  X-shell  emission,  this  is 
the  primary  type  of  radiation  produced.  The  integral  of  these  curves  over  time  gives  a  total 
yield  of  Yt  =  2.44  kJ,  Yk  =  2.01  kJ.  and  Yl  =  0.29  kJ:  this  is  a  A'-shell  efficiency  of  82%. 
The  A-shell  emission  occurs  here  only  after  pinch  as  the  plasma  expands  and  cools  down 
through  the  optimum  temperature  (~  40  eV)  for  this  radiation;  the  reason  that  an  L- shell 
pulse  does  not  occur  symmetrically  before  the  pinch  is  because  the  (artificial)  thermalization 
procedure  suddenly  jumps  the  temperature  over  this  L-shell  range.  Fig.  17(f)  shows  the  his¬ 
tory  of  the  average  core  ion  charge  state  Ze(£).  roughly  following  the  temperature  evolution. 
That  energy  is  well-conserved  throughout  the  simulation  is  evidenced  in  (g),  where  the  solid 
curve  is  the  total  energy  delivered  to  the  PRS  from  the  external  circuit  as  a  function  of  time: 


E,n{t)  =  /0  I{t')V(t')  dt'.  The  dotted  curve  Is  the  sum  of  total  kinetic  energy,  total  core 
and  skin  thermal  energies,  stored  magnetic  field  energy,  and  the  total  radiated  energy  as  a 
function  of  time.  The  final  value  of  about  5.5  kJ  is  composed  of  2.4  kJ  of  radiation  and  3.1  kJ 
of  residual  kinetic  energy  of  expansion,  the  temperatures  and  magnetic  field  having  largely 
fallen  off  (as  shown  in  (h),  where  the  magnetic  energy  \LI 2  is  plotted).  The  characteristics 
of  the  skin  are  shown  in  Figs.  I7(i-k);  the  skin  temperature  rises  to  over  80keV,  mainly  due 
to  the  fact  that  the  skin  depth  (j)  remains  so  small  ( 6  <  10~3cm)  during  the  run-in  phase, 
producing  a  large  current  density  and  hence  extreme  resistive  heating.  Actually,  the  dom¬ 
inant  heating  of  the  skin  is  due  to  resisitivity  only  during  the  first  few  nanoseconds  of  the 
run-in;  after  the  temperature  reaches  the  10-100  eV  range,  compressional  heating  (~  uT /a) 
becomes  the  dominant  effect.  Admittedly,  this  degree  of  skin  heating  is  quite  unrealistic, 
but  the  radiation  pulse  from  the  skin  (k)  is  insignificant  compared  with  that  from  the  core 
(the  skin  has  much  less  volume),  and  the  skin  only  plays  the  role  of  carrying  the  current. 
It  turns  out  in  this  case,  however,  that  the  high  pressure  in  the  skin  is  the  dominant  force 
which  bounces  the  plasma  back  out  at  pinch.  As  stated  earlier,  this  skin  effect  can  be  inves¬ 
tigated  by  varying  the  input  parameter  <§0  which  initializes  the  skin  depth.  Finally,  some  of 
the  circuit  parameters  are  shown  in  Figs.  17(1— o):  the  open  circuit  generator  voltage  trace 
(1)  is  a  piecewise-Iinear  fit  to  an  experimental  curve,  the  voltage  across  the  PRS  element 
in  (m)  reverses  sign  as  the  plasma  bounces  (due  to  the  reversal  of  the  impedance  L  —  u 
shown  in  (o)).  and  the  load  current  flowing  in  the  PRS  has  a  slight  hitch  when  pinch  occurs 
(~  35  ns  after  peak  current)  because  of  the  smallness  of  the  impedance  then  IL  ~  u  —  0). 


TEMPERATURE  k  CORE 


FIG.  17.  Typical  behavior  of  dynamical  variables  and  model  characteristics  in  optimum 
p  case,  (a)  outer  radius,  (b)  outer  radius  velocity,  (c)  density,  (d)  core  temperature. 


It  is  difficult  to  compare  the  evolution  of  the  individual  characteristics  shown  in  Fig.  17 
with  experiment  for  obvious  reasons:  the  outer  radius  in  experiment  is  not  well  defined, 
time  resolution  of  characteristics  is  not  precise,  temperature  and  density  are  difficult  to 
measure,  etc.  One  can,  however,  compare  grosser  features  such  as  circuit  quantities  and 
aspects  of  the  radiation.  In  Fig.  18  we  show  a  typical  comparison  of  (a)  load  current,  (b)  K- 
sheli  radiated  power  pulse,  and  (c)  K -shell  radiated  energy  as  a  function  of  time.  While  the 
current  appears  to  match  very  well,  the  radiation  pulse  produced  by  the  model  is  significantly 
narrower  than  the  experimentally  observed  pulse.  The  reason  for  the  difference  between  the 
widths  of  the  simulation  and  experimental  pulses  in  this  case  is  only  weakly  related  to  the 
temperature  dependence  of  the  radiation  model  used  in  the  code  (radiating  over  a  slightly 
larger  temperature  range  would  broaden  the  pulse,  but  the  density-squared  effect  would  still 
constrain  the  pulse  to  occur  only  over  10-30  ns  around  pinch  time).  The  primary  reason  for 
the  difference  in  pulse  widths  is  the  so-called  “zipper  effect” :  as  the  gas  jets  away  from  the 
nozzle,  it  tends  to  spread  radially  so  that  the  real  physical  gas  puff  resembles  more  a  hollow 
cone  than  a  hollow  cylinder.  Thus,  the  end  closest  to  the  nozzle  begins  the  implosion  with  a 
smaller  outer  radius  and  therefore  pinches  (and  radiates)  perhaps  10-20  nanoseconds  before 
the  end  of  the  plasma  farthest  from  the  nozzle  (with  larger  initial  radius).11  This  axial 
dependence  of  the  radiation  emission  broadens  the  measured  pulse,  although  as  is  shown 
in  Fig.  18(c),  the  total  energy'  produced  in  simulation  and  experiment  is  nearly  the  same 
(note  the  time  scale  only  goes  to  180ns).  The  “zipper  effect”  cannot  be  modelled  with  the 
present  code  (which  assumes  axial  uniformity),  but  a  slight  modification  which  partitions 
the  gas-puff  into  several  independent  contiguous  segments  (with  separate  properties,  such 
as  different  values  of  a0.  te.  t,  and  even  p)  is  being  investigated  in  order  to  reproduce  this 
phenomenon.12 
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E  PARAMETER  OPTIMIZATION 


Now  that  evidence  has  been  presented  which  shows  that  this  model  produces  results 
which  compare  satisfactorily  with  experimental  observations  when  simulating  current  ex* 
periments,  the  code  will  be  used  as  a  tool  to  study  parameter  regimes  not  presently  available 
to  experiment  in  order  to  investigate  optimization  of  design  and  operating  parameters.  One 
such  issue  has  been  the  effect  of  changing  the  gas  nozzle  diameter:  how  would  the  radiation 
yield  and  efficiency  depend  on  the  initial  radius  of  the  plasma?  This  is  a  time-consuming 
question  to  answer  experimentally,  as  a  new  nozzle  would  be  required  for  each  data  point: 
this  could  be  very  expensive,  if  one  does  not  have  some  idea  of  the  appropriate  sizes  to  try 
(larger  or  smaller?  and  how  much?).  Therefore,  the  following  study  was  performed:  After 
normalizing  model  parameters  so  that  code  results  and  experimental  measurements  match 
satisfactorily  as  described  above,  the  total  yield  Yt  and  the  If -shell  yield  Yk  were  computed 
(for  neon)  at  fixed  load  current  Ip  =  1.2 MA  for  different  initial  radii  a0  as  a  function  of 
mass/length  n.  The  version  of  the  code  with  skin  was  used  (in  fact,  ail  runs  were  made  with 
60  =  c/we  as  in  the  results  of  the  previous  subsection). 

The  results  of  this  investigation  are  shown  in  Fig.  19.  These  curves,  labelled  by  the  value 
of  initial  outer  radius,  clearly  show  a  trend  toward  much  higher  yields  for  smaller  values  of 
a0.  For  example,  the  A'-shell  yield  with  the  oo  =  1.25cm  nozzle  is  about  2kJ  (Fig.  19(b)), 
whereas  a  nozzle  with  ao  =  0.9  cm  is  predicted  to  produce  over  4kJ.  A  reason  for  the 
improved  behavior  for  smaller  nozzles  could  be  that  higher  densities  (smaller  radii)  at  pinch 
are  reached.  Initial  values  smaller  than  ao  =  0.8  cm  were  not  systematically  investigated, 
but  it  would  be  expected  that  the  trend  of  increasing  yield  for  decreasing  nozzle  radius 
can  not  continue:  eventually  the  run-in  phase  before  thermalization  would  be  so  short 
that  the  plasma  would  not  acquire  much  kinetic  energy,  implying  only  a  small  amount 
of  energy’  available  for  thermalization.  Indeed,  code  runs  with  very  small  ao  tend  to  end 
prematurely  because  the  outer  radius  has  become  negative:  this  “crashing*  of  the  plasma 
onto  the  axis  indicates  that  the  pressure  was  not  enough  to  bounce  the  plasma  back  out 


after  thermaiization.*  The  X-ray  efficiency  shown  in  Fig.  19(c)  peaks  at  about  80%  near 
optimum  mass  values,  but  falls  off  dramatically  at  higher  values  of  m-  This  peak  efficiency 
of  80%  is  near  to  what  has  been  observed  in  G.W1BLE  U  experiments. 1 1 

The  same  optimization  investigation  has  also  been  performed  with  the  version  of  the 
code  which  does  not  have  a  skin.  The  resuits  for  YY.  Y'k  and  ( Yk/Yt )  are  shown  in 
Fig.  20.  Again,  the  trend  toward  better  yield  at  smaller  nozzle  radius  is  apparent,  although 
the  magnitude  of  the  yields  at  different  oq  are  somewhat  higher  than  in  Fig.  18.  Also,  as 
mentioned  above,  the  optimum  value  of  fi  for  a  given  current  and  initial  radius  is  computed 
to  be  slightly  larger  in  the  version  without  skin  than  in  the  one  with  skin.  Furthermore, 
in  this  version,  almost  all  computer  runs  with  a0  <  1  cm  resulted  in  the  “crashed”  mode 
(a  <  0)  described  above.  Thus,  as  stated  above,  the  general  behavior  of  the  radiated  yields 
and  X-ray  efficiency  in  Fig.  20  are  similar  to  those  in  Fig.  18;  the  reason  for  the  discrepancies 
between  the  results  produced  by  the  two  versions  of  the  code  is  being  investigated. 

Some  of  the  results  of  this  optimization  study  are  summarized  in  the  Table  in  Fig.  21. 


*  This  condition  is  signaled  by  the  error  message  “trouble*!.”  The  remainder  of  the 
system  will  continue  to  a  normal  termination  with  the  PRS  fixed  at  its  last  non-negative 
value  of  a  and  with  u  =  0.  Actually,  there  are  many  reasons  for  such  an  occurance.  including 
a  scenario  where  the  increasing  density  as  the  radius  decreases  causes  more  radiation,  which 
cools  the  plasma,  decreasing  the  pressure,  which  allows  a  further  reduction  of  radius,  higher 
density  and  more  radiation,  etc.  It  is  believed,  however,  that  because  the  equations  evolved 
by  the  code  admit  conservation  of  energy,  this  condition  should  not  occur  and  might  be 
corrected  by  decreasing  the  timestep.  In  any  case,  the  data  produced  by  such  a  run  should 


not  be  considered. 


P  (/xg/cm) 

FIG.  19.  Results  of  nozzle  radius  optimization  study.  With  load  current  fixed  at  1.2  MA, 
opacity  —  0.25  (code  version  with  skin),  the  radiation  yield  was  computed  as  a  function  of 
mass/length  p  for  various  values  of  initial  radius  ao-  (a)  Total  radiated  energy,  Yt- 


ft  (/tg/cm) 

FIG.  20.  Results  of  nozzle  radius  optimization  study.  With  load  current  fixed  at  1.2  MA, 
opacity  =  0.12  (code  version  without  skin),  the  radiation  yield  was  computed  as  a  function 
of  mass/length  n  for  various  values  of  initial  radius  oq.  (a)  Total  radiated  energy,  IV- 


VI.  CONCLUSION 


A  simple  theoretical  model  of  an  imploding  gas-puff  plasma  radiation  source  has  been 
described,  as  well  as  its  implementation  in  the  computer  code  IMPLODE.  The  purpose  of  the 
model  is  to  be  able  to  compute  the  amount  of  radiation  produced  by  the  imploding  plasma, 
yet  it  was  intended  to  be  simple  enough  to  be  used  in  design  and  parameter  optimization 
studies.  Based  on  single  fluid  MHD  theory,  these  equations  were  averaged  over  the  volume 
of  the  gas  to  obtain  ordinary  differential  equations  for  the  outer  radius,  the  outer  radius 
velocity  and  the  average  temperature  of  the  plasma.  Results  of  sophisticated  atomic  physics 
ionization  and  radiation  calculations  were  used  to  compute  the  radiated  power  density  as 
a  function  of  plasma  temperature  and  density.  The  results  of  the  code  compare  well  with 
existing  experimental  observations,  justifying  its  use  in  predicting  optimal  design  parame¬ 
ters.  Thus,  the  main  result  is  that  a  reduction  in  gas-nozzie  radius  from  the  current  value  of 
1.25  cm  to  0.9  cm  could  perhaps  double  the  A’-shell  X-ray  yield  for  a  neon  gas-puff.  Indeed, 
this  has  been  recently  observed  in  experiment.13 
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APPENDIX  A:  CONSERVATION  OF  ENERGY 


Is  this  Appendix,  it  is  shown  that  the  total  energy  of  the  Z-pinch  is  conserved  by  the 
equations  of  motion  derived  in  Section  II,  in  that  the  total  amount  of  energy  put  into  the 
system  can  be  accounted  for.  Thus,  one  expects  the  following  to  hold 

i{K  4  4  17,  4  U/2)  4  Ve5e  +  V,R.  =  Pin  =  /$  4  I2Z.  (A.1) 

Here  again,  K  is  the  total  kinetic  energy  of  the  imploding  plasma  (governed  by  (12,16)),  Ue 
and  L\  are  the  total  thermal  energies  of  the  core  and  skin  respectively,  including  ionization 
energy  (governed  by  (25,26)),  and  VeR,.  4  V,R,  is  the  total  radiated  power  from  the  skin 
and  the  core.  The  new  quantities  here  are  the  energy  in  the  magnetic  field  ~LI 2  expressed 
in  terms  of  the  inductance  L,  and  the  power  Ptn  delivered  to  the  system  from  the  external 
circuit.  The  latter  includes  both  the  resistive  heating  I2  Z  of  the  skin  and  the  effect  of  the 
time-varying  magnetic  flux  $.  In  terms  of  the  inductance,  the  I<P  can  be  expressed  as 

J*  =  /(!/)•  =  ±L(I2y  4  LI 2  =  ( \LI2y  4  \LI2.  (A.2) 

Thus,  the  power  into  the  magnetic  field  can  be  eliminated  from  both  sides  of  (A.1),  and  the 
energy  actually  coupled  into  the  plasma  should  satisfy 


~(K  +  Ue  +  V .)  +  \\K  4  V,k,  =  \ll2  4  I2Z. 


Collecting  the  evolution  equations  (16,26)  and  using  the  definitions  (25)  for  the  thermal 
energies  in  (23),  we  have  the  following 


K  =  — r(.Vcfe  4  N.T.)  4  (t ;Jmy 
a  4  o 

Uc  =  Ia’cOc  =  — — —  .Vefe  -  veRe,  (A. 

*  a  4-  6 

u.  =  IN  A  =  -  V,R,  -r  V,r,J2. 

o  to 
-61- 

Adding  these  together,  the  pressure  (temperature)  terms  evidently  cancel  and  one  finds 
±(K  4  Ue  +  U,)  4  VeRc  *  VaR„  =  (t\/£)  -  VsVJ2.  (A. 


SOT 


•ivivS.-;- 
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Of  course,  the  more  general  equations  (16,20)  give  the  same  result  when  added.  From  the 
definition  of  the  resistance  Z  in  terms  of  the  resistivity  rj,  Z  =  frj/A.  (where  At  is  the 
cross-sectional  area  of  the  skin),  one  can  readily  identify  the  resistive  heating  term  in  (A.5) 
with  I2Z.  Therefore,  comparison  of  (A.5)  with  (A.3)  indicates  that  the  evolution  equations 
(A.4)  conserve  energy  (as  stated  above)  if  one  can  identify  the  (vJB)  term  with  the  coupled 
energy  kLI2.  That  this  is  the  case  should  be  easy  to  determine  since  the  velocity  profile, 
current  density  and  magnetic  field  are  all  known  quantities. 

Before  proceeding  to  calculate  the  (vJB)  term,  however,  we  comment  on  the  result. 


By  actual  calculation,  one  finds 


\  f  dV.v-U 


*  1)  =  HU  +  Li)!2 


Here.  Lo  is  the  inductance  of  the  plasma  annulus  inside  the  cylindrical  cavity  assuming  that 
the  plasma  is  a  perfect  conductor  (so  that  the  current  flows  in  a  skin  of  zero  thickness); 
thus,  this  is  the  usual  vacuum  coaxial  inductance  Lo  =  (2f/c2)ln(rw/a).  The  correction 
L\  to  the  inductance  represents  the  contribution  from  the  fact  that  the  skin  has  a  finite 
thickness  <5.  As  this  correction  term  is  generally  small,  however,  for  simplicity  in  the  code 
implementation  the  imploding  Z-pinch  has  been  presented  to  the  exterior  circuit  as  an 
element  with  the  coaxial  inductance  Lo-  The  (vJB)  term  in  the  equation  of  motion  (16,17) 
is  therefore  simply  set  equal  to  \LoI2  in  order  to  be  consistent.  Thus,  the  plasma  evolves 
in  such  a  way  that  the  energy  it  recieves  from  the  exterior  circuit  is  exactly  that  delivered 
to  it  as  computed  by  the  current  and  voltage  on  the  element. 

Omitting  the  details,  the  (vJB' I  term  is 

<a-j> 

__ePu  r  _b_  2 a  ab  f  2g2ln(a/j)V 
“  c2  a  a-rb  + a  +  b  a2—  z2  \  a2  -  z2  ) 

With  Lo  =  -(2C/c2)(ri/a),  comparison  of  (A.7)  with  (A.6)  gives  the  formula  for 

i,  «H[  ‘  (A.8) 

c2  a  a-rb  a  +  b  a2  -  x2  \  a2  -  z2  / 


Now,  defining  the  inductance  in  terms  of  the  energy  ~LI 2  in  the  magnetic  field,14  one 
obtains 

±I/2  s  J  dV  {B2/ 8») 

II2 ,  ,  ,  ,  2s*  r  J  n2,  , 

=  k(r»/a)  +  -g -  j  r*r  B2[r) 

=  §(Lo  +  AX)72, 

where  the  contribution  from  the  skin  gives  the  extra  inductance 


A  L 


-21  1  A 

”  C2  02  -  X2  \4 


(a2-*2)-x2  +  ^^ 


1S  —  X2  / 


The  time-derivative  of  (A.10)  is 


b  2a  ab  / 
a-ri>  a  -r  h  a2  —  x2  \ 


2x2ln(a/x)  Y 


(A.9) 


(A.10) 


(A.ll) 


c--x2  /. 

again  using  x  =  f(x).  This  is  the  same  as  (A.8),  so  that  one  can  correctly  identify  the  (vJB) 
term  as  the  inductively  coupled  energy  ~LI2. 
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APPENDIX  B:  EQUATIONS  OF  MOTION 


In  this  Appendix,  the  expressions  for  the  ordinary  differential  equations  actually  imple¬ 
mented  and  evolved  in  the  subroutine  IMPLODE  are  given.  As  derived  in  Section  II,  these  are 
five  coupled  equations  governing  the  outer  radius  a,  the  outer  radius  velocity  u,  the  average 
total  “thermal”  energy  of  the  core  Qe  and  skin  Q„,  and  the  core  radius  x  —  a  —  6. 


A.  EQUATIONS  WITH  SKIN 

By  definition,  the  first  equation  is  simply  a  =  u.  The  equation  for  the  velocity  u  now 
comes  from  the  kinetic  energy  equation  (16,17);  with  the  use  of  (7)  for  the  velocity  profile 
u(r),  the  total  kinetic  energy  K  can  be  computed  and  then  differentiated  with  respect  to 
time.  Thus, 


K  =  J  W$pv2(r) 

=  irp£j  rdr  (Afr-2  +  2AjA2  +  Ajr2)  (B-l) 

=  [A\  In (o/6)  +  A,  A,(o3  -  63)  +  \A\{d ‘  -  64)] . 

This  expression  remains  valid  even  after  the  assembly  time  (when  the  inner  radius  6  — ►  0) 
provided  the  correct  values  of  Ai  and  A?  are  used: 


A 2 


Ai 

abu  u 
a  •+■  b  o  +  6 

0 


(B.2) 


6  >  0 : 

6  =  0: 

a 

With  these  values  in  (B.l),  the  kinetic  energy  before  and  after  6  —  0  takes  the  following 
forms: 


u 

a 


b  >  0 : 
6  =  0: 


K  =  ipfu5 
K  =  \filv2. 


2ab  4a363ln(a/6) 
(a  +  6)2  +  (a  -  6)(a  +  6)3 _ 


(B.3a) 

(B.3b) 


Now,  differentiating  these  with  respect  to  time  one  obtains 


6  >  0 : 


K  -  ,Uui  [ 

A  "  2 (a +  6)2  [ 


a?  +  b2  +  Aab  + 


4  a2b2  ln(o/6) 


a3  —  b3 


fdu3 
2(o  +  6 


iu  3 

+  6)3 


+  62  -  4a6  + 


4ab(a2  +  62  - 


62 -a6)ln(a/6)] 

-  ’  <B-4a) 


6  =  0:  K  = 


(B.4b) 


where  the  constant  thickness  condition  a  =  6  s  u  has  been  used.  Note  that  as  6  — *  0,  the 
coefficient  of  u  in  (B.4a)  smoothly  goes  into  the  correct  expression  (B.4b)  for  6  =  0;  the  u3 
term,  however,  does  not  vanish  in  this  limit  but  goes  to  the  value  u3/<*.  This  discrepency  is 
of  course  due  to  the  abrupt  qualitative  change  in  the  velocity  profile  as  6  — *  0. 

The  terms  on  the  right-hand  side  of  (16)  are  now  calculated  as  follows:  the  pressure 


terms  are 


-^ThNcte  =  —nVe(Ze  +  1) 

a  +  6  cto 


=  2fUufe 


(2e  + 1)  x2  -  V2 

m,(a  +  b)  a3  —  b2’ 


7TiN‘*‘  ~  aT%nV'^'  * l) 


=  2  uluf, 

while  the  (vJB)  term  is  simply  replaced  by 


(2,  +  1)  a2  -  x2 
m<(o  +  6)  o2  —  6s  ’ 


*  **  r2 


(att)  -  ^o/2  =  -J 


(see  the  discussion  in  Appendix  A). 

Combining  these  expressions  in  the  form  (17),  the  equation  for  the  outer  radius  velocity 


u  =  -(CyCj)u2  +  (CWCi)fc  +  [CufCx)t'  -  (C4/Ci)/2, 


where  the  coefficients  are 


(B.9a) 


Ci  = 

c2  = 


a*  +  62  +  4ab 


(a +  6)2  [ 

1 


4o262  ln(o/6) 


(a +  6)3  [ 


a2  +  6'  -  4a6  -r 


a2  —  6s 
4a6(a2  +  62  -  a6)  ln(a/6) 


a2  —  6s 


c  _  4(2,  +  1)  X2  -  62 
3c  mi(a  +  6)a2-62’ 


_  4(Z«  +  1)  a2  -  x2 
3*  m^(a  +  6)  a2  -  br  ' 


C<  = 


HC?a 


While  6  =  0 


Ci  =  1, 

C2,  H  0, 

_  4(2e  +  1)  g2 


^3e  — 


mta 


^  _  4(Z,  +  l)a2-*2 

<-'3« - - - 

m,-a  a2 
/w^a 


(B.9b) 

(B.9c) 

(B.9d) 

(B.9e) 


(B.lOa) 

(B.lOb) 


(B.lOc) 


(B.lOd) 

(B.lOe) 


Again,  all  coefficients  go  smoothly  as  6  — *  0  into  their  values  at  6  =  0,  except  for  Cs-  Note 
that  these  “coefficients*  are  really  time-dependent:  they  depend  on  the  values  a,  6  =  a  —  w, 
x.  Ze(Tc)  and  Z,(T,). 

The  equations  governing  the  thermal  energies  of  the  core  and  skin  given  in  (26)  are  of 
the  form 

Qe  =  —CscVTc  +  Ctcl2  —  C7e,  (B.lla) 

Q»  =  -Cs,\lTs  +  C$,1*  —  C7,,  (B.llb) 


where  the  coefficients  are 


~  _  * 

C'4c<5'~  3(0  +  by 

(B.12a) 

P 

in 

0 

(B.12b) 

2  TO,  o2  —  62 

3»/nr(2.  +  l)  (a2-!2)2’ 

(B.12c) 

27rm^(a2  -  h2)  B 

We,7*  -  ,  ,x  "e.*- 

3m(^«,*  -r  1) 

(B.12d) 

These  expressions  are  all  valid  even  with  6  =  0.  As  prescribed  in  the  model,  no  current 
flows  in  the  core  (B.12b).  The  conductivity  e  used  in  the  skin  (B.12c)  is  taken  to  be  just 
the  Spitzer  conductivity 

f3/i 

a  =  1.96  x  9.0  x  lO13-^, 

AZf 


where  the  Coulomb  logarithm  A  is 

_  f  23.0  -  In (n2.3/2fr3/3)  for  f.  <  lOeV 

\  24.0  -  In (nl'*2y2f-1)  for  f.  >  lOeV 

These  quantities  are  computed  at  each  time-step.  The  radiated  power  density  jZc,,  for  the 
skin  and  the  core  is  computed  as  a  function  of  the  density  and  temperature  of  each  region 
(see  Section  IV). 

The  final  equation  of  motion  is  the  condition  that  the  core  radius  x  (or  boundary 
between  the  core  and  the  skin,  see  Fig.  1)  move  with  the  fluid  velocity  r(x)  at  that  radius: 


x  =  A\X~1  +  A^x, 


(B.13) 


where  the  appropriate  values  of  At  and  A?  from  (B.2)  are  used.  Combining  (B.8.B.11.B.13), 
the  five  evolution  equations  are 


u  =  -(CVCOu2  +  (CWC,)fc  +  (Cs./Ci)f.  -  (C</C,)/2, 

Qc  =  —CieUTe  +  C*/2  -  C7C, 

<5,  =  ~Cs,uT,  +  Cs»I2  -  C7,, 

x  =  A\x~ 1  4-  A-iX. 


(B.14a) 

(B.14b) 

(B.14c) 

(B.14d) 

(B.14e) 


k 


62 


These  equations  are  integrated  in  the  main  routine  IMPLODE  using  a  simple  Euler  method 
(y«+i  =  jin  A  t  4-  yn)  with  all  quantities  on  the  right-hand  sides  of  (14)  computed  at  each 
time-step  from  (B.2.B.9.B.10.B.12). 


The  current  I  flowing  in  the  plasma  radiation  source  load  element  is  computed  by  the 
transmission  line  simulation  from  the  voltage  values  on  either  end  of  the  element  divided  by 
the  load  impedance.  This  impedance  Zl  is  taken  to  be  the  impedance  of  a  purely  inductive 
eiement.  using  the  expression  for  a  coaxial  transmission  line:  Zi  —  Zq  \n{rw/a).  The  coef¬ 
ficient  Zo  for  a  standard  vacuum  transmission  line  is  (2/c),  so  that  when  the  transmission 
line  inductance  Lq  =  (2£/c~)\n(rw/a)  is  divided  by  Zi  one  obtains  the  vacuum  element 
transit  time  rL  =  (Lq/Zl)  =  (£/c).  For  implementation  in  the  transmission  line  simulation, 
however,  the  gas-puff  is  represented  to  the  remainder  of  the  circuit  as  an  element  which  is 
only  one  time-step  A t  long.  Therefore,  with  tj,  =  At  =  {Lq/Zl),  one  finds  the  appropriate 
(artificial)  value  of  Z0  to  be  Z0  =  (2/c)(£/cA0-  This  is  artificial  in  that  the  impedance  is 
larger  (or  smaller)  than  the  vacuum  impedance  by  a  factor  of  (vs/c),  where  vs  =  (£/A t) 
is  the  speed  of  a  signal  which  traverses  the  eiement  of  length  t  in  one  time-step.  For  the 
standard  GAMBLE  II  gas-puff,  l  —  4cm  and  the  usual  simulation  time-step  is  At  =  0.1ns, 
which  gives  v$  =  (4/3)c. 

Aside  from  this  technical  point  about  the  characteristic  impedance  of  the  inductive  ele¬ 
ment.  it  is  important  to  note  that  this  treatment  excludes  the  effect  of  the  energy  dissipation 
in  the  load  due  to  the  resistive  heating  of  the  skin  (B.ll.B.12c).  That  is,  the  skin  is  allowed 
to  gain  energy  due  to  I2 /a  but  this  energy  is  not  removed  from  the  circuit  by  a  correspond¬ 
ing  series  impedance  of  the  form  £  =  (£/trAs)  as  {£/2iraa 6),  where  As  =  r(o2  —  x2)  is  the 
cross-sectional  area  of  the  skin.  In  most  cases,  this  energy  sink  is  negligible  in  terms  of  total 
energy  conservation  (see  Appendix  A),  even  though  it  can  be  responsible  for  heating  the  skin 
to  anomolously  high  temperatures  (see  Section  V).  In  other  cases,  however,  the  dissipation 
is  significant  and  can  cause  large  departures  from  total  energy  conservation;  if  this  occurs, 
one  should  either  account  for  the  loss  by  incorporating  a  series  resistance  £  as  above,  or 
turn  off  the  skin  resistive  heating  term  in  (B.llb)  (t.e.,  set  Cq,  =  0  by  letting  c  —  oo).  If 
the  latter  method  is  chosen  (as  was  done  in  some  of  the  results  of  Section  V),  then  there 


is  really  no  difference  between  the  equations  (B.ll)  for  the  skin  and  core  thermal  energies. 
Thus,  one  is  led  to  consider  a  simpler  model  which  has  no  skin  as  described  below. 


B  EQUATIONS  WITHOUT  SKIN 

Since  the  magnetic  force  term  I2)  in  (B.14b)  is  computed  using  (B.7)  which  ignores 
the  presence  of  a  skin,  one  can  set  the  skin  depth  to  be  zero  implying  x  =  a  —  6  —*  a  in 
(B.9.B.10).  Now  the  equations  of  motion  are  for  just  a.  u,  and  Qe: 

a  =  a,  (B.15a) 

u  =  -(CVCOu2  -f  {C^/C^Tc  -  (C</C,)/2,  (B.15b) 

Qc  =  ~QseuTc  —  Cle,  (B.15c) 

where  the  coefficients  are  given  by 
While  b  >  0 


Ci 

C2 

&3C 

C4 


C5e 


Cl  c 


1 

a2  —  62  +  4a6 

(0-6)2 

1 

^a2  —  62  —  4a6 

(a  4-  6)3 

4(2c*l) 

m^(a  -r  6)  ’ 

2 

fic2a' 

_  4 
3 (a  t  b)  ’ 

27rm,(a2  -  62)  * 


4o262  ln(o/6) 

a2  -62  J  ’ 

4o6(a2  —  62  -  a6)ln(a/6)’ 
a2  —  62 


(B.16a) 

(B.16b) 

(B.16c) 

(B.16d) 

(B.17e) 

(B.17f) 


C  4 


While  6  =  0 


Ci  —  1, 

C3  =  0, 

„  _4(2e  +  l) 

03e  —  — — — , 


(B.18a) 

(B.18b) 

(B.18c) 


(B.18d) 


Cse~~2(a  +  b)' 

„  2irmj(a2  -  b2) 

C7C  -  3/i(fe+"l)-'Re‘ 


(B.19e) 


(B.19f) 


C.  SUMMARY  OF  VARIABLES 


Finally,  the  equations  in  Section  II  were  derived  using  Gaussian  units  even  though  exper¬ 
imentalists  measure  most  external  circuit  quantities  in  mks  units.  Therefore,  for  purposes 
of  clarification  of  several  constant  conversion  factors  that  are  to  be  found  in  the  code,  we  list 
here  the  physical  quantities  and  the  units  in  which  they  are  cast  for  purposes  of  calculation 
in  the  code  and  on  ouput  for  comparison  with  experiment: 


Quantity 


Calculation 

Units 


Output 

Units 
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